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ABSTRACT 

We describe the Cool Opacity-sampling Dynamic Extended (CDDEX) atmosphere mod- 
els of Mira variable stars, and examine in detail the physical and numerical approxima- 
tions that go in to the model creation. The CODEX atmospheric models are obtained 
by computing the temperature and the chemical and radiative states of the atmo- 
spheric layers, assuming gas pressure and velocity profiles from Mira pulsation models, 
which extend from near the H-burning shell to the outer layers of the atmosphere. Al- 
though the code uses the approximation of Local Thermodynamic Equilibrium (LTE) 
and a grey approximation in the dynamical atmosphere code, many key observable 
quantities, such as infrared diameters and low-resolution spectra, are predicted ro- 
bustly in spite of these approximations. We show that in visible light, radiation from 
Mira variables is dominated by fluorescence scattering processes, and that the LTE 
approximation likely under-predicts visible-band fluxes by a factor of two. 

Key viTords: stars: variables: Miras - stars: AGB and post-AGB 



1 INTRODUCTION 

Mira-type variable stars represent the last easily observable 
' stage in the evolution of solar-mass stars. Due to the inter- 
, actions between pulsation, shocks, complex chemistry and 
1 radiation pressure, the environment between the continuum- 
forming photosphere and the dusty wind is complex and 
difficult to model. For this reason, it has not been possible 
to extract fundamental parameters of Mira variables (e.g. 
mass, metallicity and mass-loss rate) from their spectra and 
light curves alone. Similarly, there are many easily observ- 
able properties of Mira variables (such as their visible light 
curves) that have yet to be put in a firm theoretical con- 
text. This is in contrast to Cepheid and RR Lyrae variables, 
where non-linear effects observable in their light curves have 
been used to d erive accurate (but model-dependent) masses 
and d istances l|Keller fc Woodll2006l : Im arconi fc Clementinil 
l2005l ). 

Furthermore, there has been little detailed theoretical 
work on the regions of M-type Mira atmospheres between 
the continuum forming photosphere and the radii approxi- 
matcly 10 times more distant where standard l|Draine fc 
[l984) dust types are stable. These layers are particularly 
crucial for understanding the inte raction between pulsation 
and mass loss for Mira variables l|Woodlll979l : lBowenlll988l : 
iHofner et a"l]|l998l ). 

There is a small set of models from other groups that 



include the effects of interior pulsation by introducing an ar- 
tificial piston at the base of the atmosphere, although lumi- 
nosity variations at this position are ignored. These calcula- 
tions have produced physically consistent dynamical models 
of the photosphere combined with self-consistent chemical 
equilibrium calculations and radiative transfer. For C-type 
Mira variables, this kind of literature i s most ex tensive, with 
the paper series beginning with .Hofner et al.l (|1998) now 
including a full treatment of homogeneous nucleation and 
non-grey radiative transfer as part of the 1-D dynamical at- 
mosphere code. 

For M-type Mira variables, models are faced with a 
more complex chemistry and dust formation, and at least 
for Miras with periods less than about 500 days, chaotic 
motions of the atmosphere play a larger role in dynamics 
than the winds resulting from radiative acceleration of dust 
grain s. Several app r oache s have been attempted. The mod- 
els of lHofner et al. use a non-grey radiative transfer 
in their dynamical models with mean opacities in 51 fre- 
quency meshes, do not consider the formation of dust grains 
in detail and h ave not yet been co mpared with observations. 
The models of Ijeong et all (|2003h use a grey approximation 
for molecular opacities and a composition-independent dust 
opacity, but have a sophisticated treatment of the dust nu- 
cleation processes. This approach is more applicable to the 
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longer period very dusty Mira variable they modelled than 
to the more common Miras with periods less than 500 days. 

There are no 3-D models of Mira variables at this time, 
but 2-D models of C-rich Mira variables by i Woitkc (2006) 
show the key properties of the shock-dominated dynamics: 
chaos over large spatial scales similar to the cycle-to-cycle 
variations seen in 1-D models, and fine structure in the shock 
fronts caused by the Ray leigh- Taylor instability. 

The models in t his paper are an improvement on the 
iHofmann et al.l (|l998l ) modelling scheme. Those models were 
one-dimensional and were based on a three-step modelling 
process to optimize computational efficiency. First, a grey 
dynamical model was constructed in which pulsation was 
self-excited: i.e. there was no 'piston' artificially causing the 
pulsation, and luminosity variations arising in the interior 
were accounted for. Second, a non-grey radiative transfer 
scheme was used to calculated the detailed temperature pro- 
file of the atmosphere. Finally, more detailed integration 
through the final atmospheric profile enabled spectra and 
center-to-limb profiles to be calculated. 

Here we describe our Cool Opacity-sampling Dynamic 
Extended (CODEX) atmospheric model series for M-type 
(oxygen-rich) Mira variables. The models include self- 
excited pulsation with new approximat ions for convective 
energy transport (|Keller fc WoodI |2006| ). and an opacity- 
sampling method for radiative transfer in LTE. These mod- 
els are described in detail in Section O They are tested first 
in Section |3] by further calculations that examine the errors 
in the approximations caused by our three-step modelling 
process, and some exploratory calculations into the effects 
of the LTE approximation and the effect of the dynami- 
cal atmosphere (i.e. velocity stratification) on the radiative 
transfer. The models are further tested in Section |3] by com- 
parison with spectroscopy and optical interferometry, with 
the aim of determining which model predictions are expected 
to be most reliable. 



The self-excited pulsati on models were made w ith the 
pulsation code described in iKeller fc Wo^ { 20061 ) . Given 
the physical input parameters above, and adopted values of 
the model parameters a^n (mixing-length in units of pressure 
scale heights) and (the turbulent viscosity parameter) the 
model naturally pulsates, with the amplitude limited by the 
non-linear loss processes of shocks in the outer atmosphere 
and turbulent viscosity. After beginning with a static model, 
the nonlinear pulsation model was run to limiting ampli- 
tude. The period of this model was then compared to the 
desired period of 330 days, and the amplitude was compared 
to the observed pulsation amplitude (from the photometry 
of Whitelock et al . 2000). The value of am was then ad- 
justed until the correct period was obtained, and was 
adjusted to give the correct pulsation amplitude. The final 
values adopted were am = 3.5 and a^ = 0.25. The static 
model of the so-called parent star with 5400 has a pho- 
tospheric radius Rp — 215 R© and an effective temperature 
TcS = 3380 K. The pulsation period of this model in the 
linear approximation was 330 days. 

The long-term behavior of these models is illustrated 
in Figure [T] The shaded regions are the phases chosen for 
input into the more detailed atmospheric models. In this 
paper we only consider models from the third shaded region 
(between ~7000 and 8000 days), which we label the fx series 
of models. 

The grey radiative transfer in the ou ter layers of these 
models is based on Rosseland op acities of [Rogers fc IglesiasI 
(jl992)) and lFerguson et al.l (|2005l ). In the outer atmosphere, 
the simple ad-hoc modification of th e radiative diffusion 
equation at small optical depths used bv lFox fc WoodI (|1982| ) 
is adopted. This approximation guarantees T cx r~°'® in the 
outer layers. However, this approach is not accurate enough 
for spectrum computation, so we choose to re-solve for the 
gas temperature, as described in the following section. 



2 MODELLING ASSUMPTIONS AND 
METHODS 

2.1 The pulsation models 

The aim of the pulsation models was to produce a Mira 
model with a period of close to 330 days, which matches the 
period of the local Mira variable o Ceti. T he near-infrared 
JHK photometry of lWhitelock et al.l (|2000l ). combined with 
the distance e stimate of 107 pc to o Ceti from the Hippar- 
cos parallax (k napp et"al]l2003h or the LMC (Mboi, loR-P) 
relation for a 3 30 day Mira variable (|Hughes fc Woodlll990l : 
iFeast et al.lll989ll , yields a mean luminosity of close to 5400 
Lq. This mean luminosity was adopted for the model; giving 
it the name o54. Given the luminosity, the luminosity-core 
mass relation (e.g lWood fc Z arro 198lj) defines the mass of 
the core (0.568 M©) below the inner boundary of the pul- 
sating envelope, which is set at about lO'^ K and a radius 
of 0.3 R0. A mass of 1.1 Mq was adopted, to match the 
mass estimated f rom Galactic kin ematics for a Mira of pe- 
riod ~330 days (jjura fc Kleinman n 1992 ) . A metal abun- 
dance Z = 0.02 was adopted, along with a helium abundance 
Y=0.3, which is close to the envelope helium abundance of 
a 1.1 Mq star on the AGS after it ha s undergone first and 
second dredge- up (|Bressan et al.|[l993l ). 



2.2 The atmospheric models 

With the gas pressure fixed from the dynamical mod- 
els, the gas and dust temperatures in the photosphere 
are re-solved by using an opacity-sampling method in 
LTE. This is in contrast to previous models of Mira 
atmospheres that us ed a small 72-wavelength averaged- 
opacit y mesh (e.g. Hofmann et al.l 1 19981 : llreland et al.l 



opacit y mesn (e.g. motmann et al.l llijijal : lireiang et al.l 
2004bllal: llreland fc Scholjl2006l') or 51-wavelength mesh (e.g. 



Hofner et al1l2003l ). Our continuum and Just-Overlapping- 



Li ne-Approximation (JO LA) VO opacities are the same as 
in iHofmann et al.l (|l998l ). The continuum opacity sources 
(free-free and bound-free where applicable) are H, H-, H2, 
H2-, He- and Thompson scatterin g. The H2O lines come 
from Partridge fc Schwenkel l| 19971 ) and those of TiO from 
ISchwenkel (|l998l ). Other diatomic molecular (CO, OH, CN, 
SiO, MgH) and neutral atomic lines (Na, Mg, Al, K, Ca, Ti, 
V, Cr, Mn, Fe,and Ni) come from the input to the ATLAS12 
models (|Kurucdll99i ). The lack of metal bound-free opac- 
ities means that our models are inaccurate at wavelengths 
shortwards of 450nm: wavelengths that are both unimpor- 
tant energetically in the atmospheric layers modelled by the 
opacity-sampling code and strongly affected by numerous 
other approximations such as LTE. 

A micro-turbulence of 2.8 km/s is assumed, in bet ween 
that used in the M-giant models of iPlez et all l|l992l ) and 
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Figure 1. The luminosity (top panel), effective temperature (cen- 
ter panel) and radius of selected mass zones (lower panel) as a 
function of time for the o54 pulsation model. The red dashed line 
in the bottom panel corresponds to the radius at Rosseland mean 
optical depth 2/3, and the effective temperature (oc {L / B?Y^'^^) 
refers to this radius. 

that derived by iHinkle fc Barned (| 19791 ) . Line profiles for 
molecules are taken to be Gaussian, as pressure-broadening 
is negligible and the opacity is dominated by a very large 
number of weak lines. For atomic lines, a Voigt profile is 
assumed, with only radiative damping taken into account 
with a nominal F value of 10*, typical of atom ic transitions . 
The equation of state is based on that used in iTsuiil (| 19731 ) ■ 
with 35 atoms and their first 2 positive ionization states, 60 
molecules, H-, H2- and F-. The chemical equilibrium con- 
stants for CO, MgH, SiH and Ti02 have since been m e asure d 
more accurately and differ significantly from iTsuiil lll973fl . 
so for these molecules we have adopted constants from 
ISharp fc HuebneJ ('l990Y We as sume solar abundance s, and 
take the solar abundances from lGrevesse et al.l l| 19961 ). 

The dust opa city and equation of stat e is taken from the 
approximations of llreland fc Schol3 (|2006l ) . with parameters 
designed to best fit dust scattering opacity observations: the 
log of the number of nuclei per H atom log(7Vnuc) = —12.8 
and the ratio of sticking coefficients as = 1. This treatment 
of dust opacity is temperature-dependent, due to the ~ l/xm 
dust opacity being a strong function of the Fe-content of the 
dust, only approaching stan dard interstellar medium dust 
opacities (e.g. lOraine fc Leelll984l ) at T < IQOOK. 

The opacities are calculated for 4300 wavelengths run- 
ning from 200 nm to 50 /im, with 1 nm sampling between 
200 nm and 3 /xm and sampling proportional to for longer 
wavelengths. The effect of the velocity stratification on the 
radiative transfer is not taken into account, but is expected 
to have significant effects only in the vicinity of shock fronts 
(see Section |3]). 

Instantaneous relaxation of shock heating behind 



the moving shoc k front is assumed (see discussion in 
iBesseU et alll989l ). This is equivalent to enforcing L — Lsurf 
everywhere, with Lsurf the model surface luminosity. Al- 
though the shock luminosity is important for modelling 
emission lines, its luminosity is generally negligible in the 
photospheric regions where the continuum is optically thin 
and where lines are formed. At phases immediately before 
maximum luminosity, as the light curve of a Mira is rapidly 
rising, the shock luminosity can for a small time reach half 
the model luminosity. The assumption L = Lsurf is incor- 
rect in this case, meaning that the shape of the continuum 
during this rapid pre-maximum flux increase phase may not 
be a good approximation. The effects of this approximation 
are examined further in Section [3.31 

The equa tion of radiative transfer is sol ved using the 
same code as ISchmid-Burgk fc Schofl { 19841 ). with up to 
80 discrete depth points. Since this code uses a spline in- 
terpolation, the pressure discontinuity at the shock front is 
smoothed out over 0.05-Rp (deep layers) to 0.1-Rp (high lay- 
ers). The model outer-boundary is fixed at 5Rp: a boundary 
where the dust formation is more complex than our simple 
approximations allow, and where radiative acceleration on 
dust can influence the dynamics signiflcantly. Layers outside 
5 Rp can be considered the 'wind' zone, and can be seen ob- 
servationally only in the mid-infrared due to dust opacity 
and in the strongest of molecular transitions. 

Once the equation of radiative transfer is solved, the 
spectrum and center-to- limb proflle is computed over a flner 
wavelength grid, typically with 17000 wavelengths. This 
flner grid includes wavelength in the radio regime, where 
the wavelength-dependence of free-free opacities means that 
the continuum phot ospheric radius i s quite different to that 
in the near- infrared iReid fc MentenI l|2007h . 



3 DISCUSSION OF MODEL QUALITY 

3.1 Verification of opacity-sampling code for a 
static model 

Previously published models for sta tic extended M- 
giant atmospheres fro m other groups (|Fluks et al.l Il994l : 
iHauschildt et al.lll999l ) have demonstrated relatively good 
agreement with observations. Therefore, as a flrst step in 
verifying the performance of our code, we will compare the 
model spectra for a relatively compact static model to those 
from the PHOENIX c ode. This code has also adopted LTE for 
published models (|Hauschildt et al"] Il999l ) , so if the input 
opacities to their models and our models were the same, 
one would expect the model outputs to be the same. As 
a comparison model, we choose a 5 solar-mass model with 
log((/)=0 (cgs units), and Toff = 3200 K. This mode l is one 
of the most extended conflgurations from Hauschil dt et al.l 
(| 19991 ). which is not designed for very extended atmospheres, 
and is the most compact model where our numerical approx- 
imations still solve the equation of radiation transport with 
reasonable accuracy. 

Figure [2] shows the model spectra comparison between 
1 and 3.0 /xm, and Figure |3] shows the comparison between 
0.5 and 1 f^m. The most noticeable difference is the depth of 
the TiO features: the dominant molecular bands short- wards 
of 1 ^m. Some of this difference comes from the use of the 
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|j0rgensenl l|l994 ) line list by iHauschildt et all (|l999l ) (N.B. 
more recent PHOENIX model s of M-dwarfs us e the same TiO 
opacity as we do) . When the |j(zirgensenl (|l994h line list is used 
with the CODEX models (adopting the oscillator strengths 
from lAllard et al1l2000l ). the difference at around 750 nm 
becomes smaller and the details of the band shapes agree 
more, as seen in Figure|3] The discrepancy in the TiO feature 
at 1.25 /im as shown in Figure[2]is also reduced using this line 
list (not shown). However, the overall flux level of models the 
at shorter wavelengths, particularly between 500 and 700 nm 
remains discrepant. Both reducing the effective tempera ture 
of our models by ~100 K and using the |j0rgensenl (| 19941 ) line 
list produces a very small difference between the CODEX and 
PHOENIX models (not shown). 

We also compare the 3200 K static model to observa- 
tions in Figure [5] where an M6 and an M7 spectrum is 
compared to the models. We note that the difficulties in 
our models fitting the gap between the strong TiO bands at 
750 nm is in com mon with the curren t version of the PHOENIX 
models (e.g. see lLancon et al]|2007l ). due to an inadequacy 
in the line lists. Until this line list is updated, spectral in- 
dices based on TiO absorption will not be able to be used 
to convert the model spectra to spectral types. 

At ~ 2.5^m, the PHOENIX models have a feature that 
is not in our models: presumably this is due to the differ- 
ence in water line lists: again the line list used by us is that 
adopted in more modern PHOENIX dwarf models and not the 
published PHOENIX giant models. There is also a difference in 
the temperature structure of the models, as seen in Figure S) 
We have independently checked the temperature structure 
of our models using a (computationally slow) Interactive 
Data Language (IDL) radiative transfer code that uses the a 
modified Unsold- Lucy temperature corre ction for spherical 
atmospheres (e.g. IHauschildt et all [l99^ and trapezoidal- 
rule integration, finding differences of order 10 K only, so we 
therefore expect the difference in temperature structure to 
be due to further opacity differences rather than computa- 
tional errors. 

Not e that previou s versions of our code, such as that 
used in llreland et al.l (2004b, a) could not reproduce the 
spectra of these static models nearly so well, due to the 
inaccuracies of the Just Overlapping Line Appr oximation at 
these low to rn oderate TiO optical depths (e.g. lBrett|[l990l : 
iTei et al.|[2003h . 

3.2 Critical assumptions examined at 3 model 
phases 

There are a large number of computational checks one can 
make on the quality of our atmospheric models before com- 
paring them directly to observations. The checks that we 
consider most relevant will be presented here, relating di- 
rectly to model physical assumptions and approximations: 
1) Is a grey temperature-profile adequate for describing the 
atmospheric dynamics?; 2) Does the velocity stratification 
significantly influence the solution to the radiative transfer 
equation?; 3) How much does the LTE approximation affect 
the spectrum and temperature profile?; and 4) How does 
neglecting the shock luminosity by insisting "L = Lsuri" ev- 
erywhere affect the model properties? We will attempt to 
answer these questions by examining three model phases in 
detail. The 3 models have phases of 0.20 0.59 and 0.80 rela- 
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Figure 2. The spectrum between 1 and 3.0 fim produced from 
the CODEX atmosphere models at 3200 K (solid lines), compared 
with that from the PHOENIX atmosphere code (red dotted line). 
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Figure 3. The spectrum between 0.5 and 1.0 fim produced from 
the CODEX atmosphere models at 3200 K (solid lines), compared 
with that from the PHOENIX atmosphere code (dotted line), and 
the CODEX models using the same TiO line list as PHOENIX (red 
dashed line) (see text). 



five to (estimated) maximum bolometric luminosity, which 
is assigned phase 0.0. 

Using a grey model for the non-linear dynamical mod- 
els is only valid to the extent that the temperature pro- 
file is consistent with that from a more accurate non-grey 
model. Furthermore, the grey model uses approximations for 
the spherically-extended atmosphere that also affect the way 
that radiative acceleration is calculated. Figure [6] shows the 
grey and non-grey temperature profiles, and the grey and 
non-grey accelerations as a function of radius. The total ac- 
celeration includes a gravitational term, a gas pressure term 
and a radiative acceleration term: 



. 1 dPga: 

a = g+ 
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Figure 4. The temperature profile of the CODEX models for a 
3200 K model (solid line) plotted with the the PHOENIX models 
(dashed line) for the same model parameters. 
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Figure 5. The spectrum between 0.5 and 1.0/jm produced from 
the CODEX atmosphere models at 3200 K (solid lines), compared 
with an M6 giant spectrum (upper d otted line) a nd M 7 giant 
spectrum (lower red dashed line) from lpTuks et ahl (Il994l) . 



Here is the extinction coefficient at frequency u, Fi, is 
the monochromatic flux per unit frequency and other terms 
have standard meanings. 

This figure demonstrates first that the grey approxi- 
mation only affects temperatures in the outer layers at the 
10-20% level (excluding shock fronts). A 10% increase in 
the grey model temperature would cause a 10% increase 
in the pressure gradient for the same mass stratification. 
This in turn changes the acceleration (both continuous ac- 
celeration and the (5-function acceleration at a shock-front) 
at the 10% level. During a pulsation cycle, the bulk of the 
non-gravitational acceleration of a given layer occurs close 
to the continuum-forming photosphere, where the diffusion 
approximation as used in the grey model is most accurate. 
The gas in the outer layers follows near-ballistic trajectories. 



Therefore, we estimate that geometric pulsation amplitudes 
for layers between 1 and 4 7?p are accurate to <10%. 

As a dynamic atmosphere has a velocity stratification 
v{r), the radiation field seen from a test point at a spectral- 
line wavelength may strongly depend, as a consequence of 
Doppler shifting, on the depth of the layer of origin and on 
the direction of the line-of-sight. There are two classes of 
v{r) effects: a geometric v{r) projection effect that is often 
moderate except for the very most extended atmospheres; 
and more significant effects at and near a shock front where 
v{r) is discontinuous at and has a strong gradient near the 
position of the shock. Sample wavelengths that, e.g., are lo- 
cated at line frequencies within a forest of blanketing lines 
below the shock appear in the continuum above the shock, 
and vice versa. This may lead to substantial errors in com- 
puted equilibrium temperature when an opacity sampling 
technique is used for treating radiative transport that is 
dominated by line blanketing. Solutions for the case of a ve- 
locity stratification v(r) have been developed for the opacity 
distribution technique, though for only rather simple cases 
of v( r) (e.g. iBaschek et al.ll200ll : IWehrse et "alll2003l : [Castoil 
|2004| ). but no solutions are known for the opacity sampling 
method (which necessarily has large gaps in-between indi- 
vidual wavelengths that are sampled). 

We have, however, numerically examined the equation 
of radiative transfer in the dynamical model atmosphere 
without attempting a full solution. This was accomplished 
by calculating the luminosity at each layer of a model strati- 
fication on a grid of ~ 10^ wavelengths, self-consistently tak- 
ing into account the velocity of each layer, by red-shifting or 
blue-shifting the line profiles according to the difference in 
projected velocities. This code was far too slow to attempt 
a solution to the equation of radiative transfer, but could 
examine the departure from the L — L^uii criterion. Fig- 
ure [7] shows a maximum 2 % deviation from the L = Lsurf 
criterion, demonstrating that the effects of the dynamical 
stratification are a relatively minor effect in calculating the 
temperature structure. The effects of the dynamical struc- 
ture on the overall spectral shape was also examined and 
also found to be minor, as shown in Figure [S] 

Deviations from LTE tend to increase with decreasing 
densities and pressures. With typically very low pressures 
(< 1 dyncm^^) in the atmospheric regions of Mira variable 
stars where spectral features are formed, these stars are a 
long way from LTE. This is because collisions can not ther- 
malise the level populations of atoms and molecules. Non- 
LTE effects are most noticeable for electronic transitions 
of atoms and TiO, where typical Einstein A coefficients of 
~ 10^s~^ greatly exceed coUisional rates of order 10^ s~^. 
Note that this is quite different to the CO molecule, where 
Einstein A coefficients for vib rational transitions are of order 
10^ s"^ l|Chandraet al.lll996l ). 

These non-LTE effects for TiO are much more impor- 
tant in Miras than in static M giants, despite static M giants 
also having Einstein A coefficients much higher than coUi- 
sional rates. This is because static M giant atmospheres are 
not nearly as extended as Miras, so the thin outer layers 
are much warmer for the same continuum-forming temper- 
ature. In turn, this means that the band depths calculated 
in LTE are much deeper for Miras, increasing the potential 
for non-LTE effects to be important in modelling. 

A complete band-model non-LTE treatment of TiO 
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Figure 6. A comparison between the grey dynamical model temperature profile and total acceleration as a fraction of gravity (red 
dashed-lines) with that from the non-grey opacity-sampling code (solid lines), for the fx series models at phases 0.20 (left), 0.59 (center) 
and 0.80 (right) from the same cycle. The 'spikes' in the acceleration are near the location of shock fronts and are an artifact of smoothing 
the shock front in the dynamical grey and non-grey models. The dominant shock fronts are mar ked with vertical dott ed lines, and the 
dotted lines in the top panel are temperatures calculated with the 72 mesh wavelength grid as in [Ireland et all (|2004d l. 



opacity would require a complete new radiative transfer code 
(both codes used in this paper tabulate opacities as a func- 
tion of LTE temperature and pressure). However, we have 
attempted to treat the TiO absorption with a partial band- 
model non-LTE treatment in order to examine some effects 
of the LTE approximation. 

We make an assumption that the electronic singlet and 
triplet ground states for TiO are populated in thermal equi- 
librium (i.e. in both rotational and vibrational equilibrium), 
and that transitions to and from excited electronic states 
occur via scattering processes. Both resonance and fluores- 
cence scattering processes are included. This assumption 
is reasonable because radiative transitions between vibra- 
tional levels are dominated by direct ro-vibrational transi- 
tions for radiation field temperatures of ~1500 K (applicable 
to the region where the strong TiO features form) . These ro- 
vibrational transitions see a thermalised radiation field due 
to strong mid-infrared transitions of other molecules (no- 
tably H.2O). Radiative transitions between electronic excited 
states (multi-photon absorption) is not as important here as 
it is for many atomic states because of the relatively large 
gap between the ground and first electronic excited state (i.e. 
in layers where the strong features are formed, the molecule 
spends the vast majority of its time in the electronic ground 
state). Collisions with other atoms and molecules are also 
negligible in the regions where the strongest features are 
formed. 

The scattering process is added to the equation for the 
source function as follows: 
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Figure 7. A comparison between the luminosity as a function of 
radius between the opacity-sampling code (straight dashed lines) 
and full dynamical modelling taking into account all wavelengths 
and the velocity stratification (solid lines). The models, from top 
to bottom, are at phases 0.20 (7623 L©), 0.80 (4092 Lq), and 0.59 
(2548 L0). 



S, B, J and k have their usual meanings of source func- 
tion, Planck function, mean intensity and absorption coeffl- 
cient (cross-section per unit mass). The isotropic scattering 
matrix cr(A, A') describes scattering from wavelength A' to 
wavelength A, in units of cm'^/g/nm. For resonant (i.e. non- 
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fluorescent) scattering, a(A', A) — 5{\ — A')cr(A), with S the 
Dirac (5-function. 

The results of this calculation is shown in Figure |8l It is 
clear that fluorescence scattering in general dominates the 
radiative processes in the l/-band, and tha t interpretations 
such a s the simplified LTE interpretation of lReid fc MenterJ 
l| 19971 ) offer a gross over-simplification of the physical pro- 
cesses that create the visible light-curves of Mira variables. 
This statement is true both in the strong absorption bands 
and in the regions of weaker absorption in-between the 
bands. Note that a similar fiuorescence approach would have 
to be applied to VO in order to model the 0.75 to 1.2 micron 
region more accurately. 

The models still do not well- represent the deep TiO 
features, as will be discussed in Section A possible expla- 
nation for this discrepancy is inaccurate data for Ti02, that 
could enable TiO to be removed fro m the gas at higher tem - 
peratures (e.g. Ilreland et al.l 120051 : [Ireland fc Schold 120061 ) . 
or non-equilibrium processes in the chemical reactions of 
TiO. An example of such a non-equilibrium process is the 
lowering of the vibrational temperature of TiO as the mid- 
infrared regions of the spectrum become optically-thin and 
the justification for the fluorescence scattering approxima- 
tion above no longer work. This would drive the reaction 
TiO-|-H20^Ti02-l-Il2 to the right, removing TiO from the 
gas. This kind of complexity is beyond the capabilities of 
the code developed for this paper. 

3.3 Other model quality issues 

The opacity-sampling method for radiative transfer calcula- 
tions is only valid if a large enough number of wavelengths 
are used. For several representative test models with differ- 
ent stellar parameters (in particular different effective tem- 
perature), we have solved the radiative transfer equation 
based on the here adopted 4300-wavelength mesh, as well as 
meshes with 2152 and 8606 wavelength samples. Differences 
of temperatures never exceeded a few tens of degrees in high 
layers for the 4300 vs. 2152 case, and were negligible for the 
4300 vs. 8606 case. Typically, including dust absorption that 
makes the atmospheric high- layer opacity "greyer", dimin- 
ishes 4300 vs. 2152 differences. We also increased strongly 
heavy-element abundances, i.e. the effects of line blanketing 
upon the temperature stratification, and found increasing 
and significant 4300 vs. 2152 differences whereas no 4300 vs. 
8606 differences showed up. Thus, a 4300-wavelength mesh 
is well suited to cover even quite unfavorable cases of M- 
type Mira atmospheres. Temperature errors of the order of 
a few tens of degrees are smaller than other errors such as 
those caused by the smoothing of shock fronts, the dynamic 
stratification, or no n-LTE effects. We note tha t the opacity- 
sampling models of lHofner fc AndersenI (|2007| ) with only 64 
wavelength points are not expected to have sufficient wave- 
length sampling to accurately reproduce the temperature- 
profile. However, those models also used a 64-wavelength 
opacity-sampling method for the dynamic atmosphere com- 
putation, which may be preferable to our grey dynamical 
models. 

As seen in Figure |6] the temperature profiles of the 
CODEX models differ noticeably fr o m the profiles computed 
from the models of Ilreland et all l|2004bl ). where the opac- 
ities were mainly approximated by mean opacities over 72 



wavelength bins and dust formation was not included. We 
see the current models as a major improvement over those 
models. This older style of opacity computation was inter- 
mediate between the current opacity sampling method and 
a grey atmosphere, as can be seen from the temperature 
profile over the range 1000 to 2000 K. At the lowest temper- 
atures (< 1000 K), the CODEX models are warmed more than 
these older models because of the presence of dust. 

In addition to limitations evident from the tests above, 
the CODEX models have clear limits where unknown details of 
complex heterogeneous dust formation become important, 
and could not be used to model e.g. OH /IR stars. These 
approximations are discussed in detail in Ilreland fc Schol3 
(2006). Importantly, dust is assumed to form in equilibrium, 
which is only applicable to the hottest (>1000K) dust. 

Finally, we have performed test calculations where the 
L = Lsurf approximation is removed, and the shock lumi- 
nosity is added to the models over a smoothed region of 
size Ar/r --^ 5%. The only test model with a noticeable 
(i.e. above numerical errors) change in model structure was 
the phase 0.80 model, which had a 540 Lq shock within the 
continuum-forming photosphere and a total luminosity of 
4090 Lq . The model that included the shock luminosity had 
5% less H-band fiux and 10% more V-band fiux. As this 
phase is at the time where the model luminosity (and light 
curve of observed Miras) is rapidly increasing prior to max- 
imum luminosity, these differences are not very significant. 



4 RELIABLE OBSERVATIONAL 
PREDICTIONS 

Given that Section |3] demonstrated that there are some ap- 
proximations that are not physically realistic, we will exam- 
ine which model outputs are most reliable for comparison 
to observations. Clearly, only those regions of the spectrum 
that are well-approximated by LTE should give reasonable 
agreement to observations. The phase, cycle and parameter 
dependence of predictions will be examined in a forthcoming 
paper II. 



4.1 Spectra 

Figure 1^ shows some observed spectra compared with mod- 
els. Instead of spectra of o Ceti, that were not readily avail- 
able in electronic form, we chose to compare the models 
with R Cha, a Mira with a 335 day period: almost exactly 
the same as o Ceti. For phases 0.2 and 0.59, an observed 
spectrum is available at a close enough model phase to give 
reasonable spectral agreement between the CODEX models 
and observations. In particular, the near-minimum spectrum 
is much better represented by the CODEX model than by 
the pr evious generation of models as examined in lTei et al.l 
(|2003l ). The infrared features of H2O and CO are repro- 
duced reasonably, but the TiO features remain too deep in 
the CODEX models. Although some of this discrepancy is ex- 
plained by the effect of the LTE approximation, the test 
calculation results shown in Figure |8] demonstrate that a 
fluorescence scattering approximation does not noticeably 
affect band depths long-wards of 750 nm. Therefore, we can 
not consider the band depths of strong TiO features to be a 
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Figure 8. A comparison between the opacity-sampling code spectrum (solid lines) between 500 and 950 nm and that resulting from a 
fluorescence scattering non-LTE treatment of the opacity (dotted lines). The dynamical atmosphere calculation is also shown (red dashed 
lines), demonstrating that this effect is very small compared with non-LTE effects. 



reliable prediction, with or without the addition of fluores- 
cence scattering. The infrared features, particularly those 
short-wards of 2500 nm, are expected to be a reasonably 
good model prediction. 

The visible-light spectrum and flux is too small to be 
seen on the linear scale in Figure |5] As visible light curves 
are one of the richest data sets for Mira variables, theoretical 
predictions for these curves as a function of physical param- 
eters is a major goal of our effort. As visible wavelengths 
are in the Wien part of the Planck function at temperatures 
relevant to Miras, fluxes and spectra are very sensitive to 
changes in the atmospheric stratification. The effects of the 
LTE approximation as discussed in Section [3 . 21 cause visible 
fluxes to be underestimated by a factor of ~2. The difficulty 
in accurately treating dust can give another few tens of % 
error. So we expect the model-predicted visible light curve 
to have uncertainties of approximately 1 magnitude. 



4.2 Diameters 

In Figure 1101 we show the range of model-predicted diam- 
eters for o Ceti for the three phases examined. In order to 
make this plot, we calculated visibility curves for filters with 
1% fractional bandwidth and calculated the best fit uniform- 
disk visibility curves aX V = 0.6 and V = 0.2. The fltted 
angular dia meter can be a strong function of where this fit 
is made (e.g Ireland et ahl 2004lJ) . We assume a distance of 
100 pc after Ivan Leeuw en (2007') for this plot. W e also show 
mea sured diameters from Irel and et al.l (|2004d ) (< 1 /xm) 
and IWoodruff et al ] ([ioOSl) (1-4 Mm) We do not show mid- 
infrared diameters from Weiner et al.l lj2003| ). because of the 
complexity of accounting for the over-resolved dust emission 
from the wind. The error bars represent both observational 
error and the differences amongst observations at different 
phcises. The J-band diameter range of the models match the 
observations very well, showing that the continuum diame- 
ters are in good agreement: certainly within 10% in diame- 
ter, corresponding to 5% in effective temperature. However, 



the measured diameters at H, K and L bands are always 
close to the largest model diameters. We suspect this not 
to be an error in modelling, but instead a physical effect 
demonstrating that o Ceti is usually surrounded by more 
extended molecular layers than the present models produce. 
Changing e.g. the assumed model mass may help resolve this 

small discrepancy. 

In [Ireland fc Schol j (|2006l ). we showed that the diam- 
eters of Mira variables at wavelengths shorter than 1 ptm 
were in general well-described by the model of dust forma- 
tion that we use here, consistent with Figure 1101 However, 
this statement was not true for the strong TiO absorption 
bands, in particular the band at 712 nm. The inclusion of the 
fluorescence scattering approximation as a test calculation 
in this paper greatly increases the model diameters in the 
TiO absorption bands (as seen by the dotted line) , meaning 
that, where necessary, more accurate predictions for model 
diameters short-wards of 1 /xm can now be produced. 



5 CONCLUSIONS AND FUTURE WORK 

We have presented the Cool Opacity-sampling Dynamic 
Extended (CODEX) modelling method: an improved scheme 
for constructing stellar atmosphere models that is tailored 
to Mira variable star atmospheres. This scheme includes 
self-excited pulsation and a 4300-wavelength grid opacity- 
sampling code to solve for the equation of radiative transfer. 

The models stand up to a variety of numerical tests, 
including adequate treatment of shock fronts and dynami- 
cal effects, and sensitivity at less than the 10% to the grey 
approximation used in the dynamical models. The models 
notably have a ~100K maximum difference in temperature 
profile when benchmarked against static PHOENIX models, 
which give the CODEX models deeper absorption shortwards 
of 950 nm. Test calculations using a fluorescence scatter- 
ing approximation for non-LTE TiO effects demonstrated 
a modest difference in overall energetics, but a factor of ~2 
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Figure 9. Smoothed model spectra (black sol i d line s), with mea- 
sured spectra of R Cha from lLancon &: WoodI 1 I2OOCII) over-plotted 
(red dashed lines). The plots are labeled by their model phases. 
The estimated visual phases of the observations are 0.25, 0.5 and 
0.76. 
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Figure 10. Model diameters between 0.6 and 4fim for the CODEX 
models (dashed lines) and the models with the fluorescence scat- 
tering approximation for TiO (dotted line). The range shown 
corresponds to the maximum and minimum diameters of the 
three models using two fitting techniques (see text). The hori- 
zontal axis is linear in wave-number so that both the full range 
and the TiO features can be shown. P oints with error bars show 
the range of measured diameters from [Ireland et al. I ll2004d) and 
IWoodrufI et~all ll2008l) . 



increase in V-band flux and much larger diameters in TiO 
absorption bands. 

Model predictions from the series presented here, as well 
as future work, will be be made available onlin^U- Work in 
progress includes the detailed analysis of the phase and cy- 
cle dependence of model properties, models with different 



http: / /www. physics. usyd.edu.au/~mireland/codex/ 



fundamental parameters including stars with longer periods 
and modified element abundances (metallicity, S-type C/O 
ratio, modified N abundance) , as well as further comparison 
of typical model predictions with observed features. 
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